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Abstract. We consider the dynamics of dark matter solitons moving through non- 
uniform cigar-shaped Bose-Einstein condensates described by the mean field Gross- 
Pitaevskii equation with generalized nonlincaritics, in the case when the condition for 
the modulation stability of the Bose-Einstein condensate is fulfilled. The analytical 
expression for the frequency of the oscillations of a deep dark soliton is derived for 
nonlincaritics which are arbitrary functions of the density, while specific results arc 
discussed for the physically relevant case of a cubic-quintic nonlincarity modeling two- 
and three-body interactions, respectively. In contrast to the cubic Gross-Pitaevskii 
equation for which the frequencies of the oscillations are known to be independent of 
background density and interaction strengths, we find that in the presence of a cubic- 
quintic nonlincarity an explicit dependence of the oscillations frequency on the above 
quantities appears. This dependence gives rise to the possibility of measuring these 
quantities directly from the dark soliton dynamics, or to manage the oscillation via 
the changes of the scattering lengths by means of Feshbach resonance. A comparison 
between analytical results and direct numerical simulations of the cubic-quintic Gross- 
Pitaevskii equation shows good agreement which confirms the validity of our approach. 
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1. Introduction 

Since the first experimental realization of Bose-Einstein condensates (BECs) of diluted 
atomic gases, the dynamics of matter waves in different trapping potentials has attracted 
the attention of many researchers and presently is a field of growing interest (TJ [2] . The 
dispersive properties of BEC induced by the quantum properties of a slow atomic motion, 
on one side, and the intrinsic nonlinearity of the system induced by the interatomic 
interactions, on the other side, make BECs ideal systems for exploring competing effects 
of dispersion and nonlinearity. This results in a variety of spectacular phenomena 
ranging from the existences of matter wave vortices [3] and solitons pEJ [5J E] to dispersive 
shocks [U[8]. 

In the standard mean field approximation [1] the BEC dynamics is governed by the 
Gross-Pitaevskii (GP) equation 

<9?/; ti 2 

for the condensate "wave function" ip, where ^(r) is the potential of the external forces 
trapping the condensate. In the standard GP approximation 

/(M 2 ) = <#l 2 , (2) 

where g = 47rh 2 a s /m is the nonlinear interaction parameter expressed in terms of the 
s-wave scattering length a s of atoms with mass m. In this approximation, only two-body 
collisions are taken into account. The scattering length a s , however, can be varied and 
even made equal to zero by means of Feshbach resonances. Then contribution of multi- 
particle (three-body) collisions becomes comparable with the two-body contribution and 
a more complicated form of the nonlinear term, such as 

/(h/f) = <#l 2 + /#l 4 (3) 

arises in the GP equation. Although the three-dimensional GP equation (OQ) is too 
complicated for analytical studies, one can often consider experimental low dimensional 
settings, e.g. cigar-shaped traps with a strong radial confinement, for which the radial 
motion of the BEC is "frozen" and Eq. ([1]) reduces to an effective ID cubic-quintic 
generalized nonlinear Schrodinger (NLS) equation [9], 

#; = -^ + ^W + /(H 2 )^ (4) 

(here standard non-dimensional variables have been introduced, holding the previous 
notation for convenience). In this reduction, the function f(p),p = \tp\ 2 , preserves 
its form Q up to a renormalization of the constants g and h. Also note that the 
conservation of the norm of the condensate wave function (normalized number of atoms) 
is guaranteed for arbitrary functions f{p)- More complicated form of the function f(p) 
can arise if we take into account the radial motion of BEC (see, e.g., [TQl [TTJ) . Thus, we 
arrive at the problem of ID dynamics of BEC with generalized form of the nonlinearity 
f(p) and confined in the axial direction x by the trap V(x). 
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As is well known, in case of f(p) given by Eq. ([2]) with g > (repulsive two body 
interaction) and under the assumption of negligible changes of the potential V(x) in the 
region of x under consideration, the NLS equation (SJ) , (J2J) , has a dark soliton solution 
corresponding to the motion of a hole in the density distribution function, p(x — vt), 
along a uniform distribution p = po with constant velocity v which depends on the 
depth of the soliton [12]. If the potential V(x) is a slowly varying function for distances 
of the order of the soliton's width, then the coordinate of the soliton is well defined and 
we can speak about value of the potential V(x) at the "point" x of the soliton position. 
In this case we can formulate the problem of motion of a dark soliton "in the potential" 
V(x) in terms of its collective coordinate. This problem has been addressed in several 
papers (see, e.g., [131 El ESI ESI EH EH EH]) where it was shown that the dynamics of 
dark solitons is quite nontrivial. In particular, if a BEC described by the standard GP 
equation (T5]) is confined in a harmonic axial trap, 

V{x) = l -to x\ (5) 
then the dark soliton oscillates inside the condensate with the frequency [13] 

This result is quite surprising because apparently it contradicts to a simple consequence 
of the Ehrenfest theorem which says that the center of mass of the whole BEC always 
oscillates with the trap frequency ujq. The remarkable theoretical prediction ([6]) was 
confirmed in the experiment [20]. Equation (El) then implies that the motion of the 
dark soliton must be accompanied by a deformation of the density distribution which 
re-tunes the oscillation of the whole condensate to the trap frequency. In this situation 
the motion of the dark soliton becomes "decoupled" from the motion of the center of 
mass and the Ehrenfest theorem is not violated. 

However, the result is not universal. If the nonlinearity function is given by 

Hp) = hp 2 (7) 

with h > 0, then the frequency u of oscillations of a shallow dark soliton coincides with 
the trap frequency p3], 

to = u (8) 

and for deep solitons this frequency is given by a more complicated expression (see [16] ) 

f /3 1 ~ V2 

uj = I 1 + — - — = } uo « 0.6572w . (9) 

This observation poses the problem of the dark soliton motion in the case of a more 
general form of the nonlinearity. 

The aim of the present paper is to address this problem by considering a non- 
uniform cigar-shaped BEC with nonlinear interactions modeled by an arbitrary function 
f(p) of the density. In particular we consider the physically relevant example of the 
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cubic-quintic nonlinearity function in ([3]), modeling the mean field effect of the two- 
body and three-body collisions in the GP equation. In this case we find that the 
relative contribution of the two nonlinear terms to the soliton oscillation depends on 
the background density p : if a deep soliton oscillates in vicinity of the maximum p of 
the density distribution, then for \h\po <C g its frequency must be equal to (jSJ), and for 
\h\Po ^ 9 to (Q. The transition from one limiting case to another is of considerable 
theoretical interest. We derive an analytical expression for the oscillation frequency of 
a deep dark solitons which reveals the explicit dependence of the oscillations frequency 
both on the background density and on the strength of the nonlinearity parameters. 
Our formulae reproduce the results considered in [15], [16] for the particular cases (121) 
and (JTj) of the nonlinearity for which the above dependence disappears. 

We also consider the case of a quintic attractive nonlinearity (h < 0) with an 
overall repulsive interaction. In this case for the density p > g/(2\h\) we get modulation 
instability of BEC, i.e. the sound velocity 

c s = vW'(Po) (10) 

becomes imaginary so that the notion of "dark solitons" looses its sense. We show 
that in this case the oscillation frequency w as a function of po becomes singular at 
Poc = fl'/(2|/i|). The form of this singularity may provide further advances in the 
understanding of the dynamics of dark solitons in BEC. 

Our theory will be based on the elegant approach of Refs. [151 EB] where the dark 
soliton was considered as a "Landau quasi-particle" whose motion in the stationary 
potential V(x) is governed by the quasi-particle's "energy conservation" law in close 
analogy with the Newtonian particle motion. The generalization of this approach to 
arbitrary f(p) allows us to derive the equation of motion of the dark soliton for the 
general case and to show the dependence of the soliton effective mass on the system 
parameters. A comparison between analytical results and direct numerical simulations 
of the cubic-quintic GP equation shows a very good agrement which further confirms 
the validity of our analysis. 

We finally remark that the dependence of the dark soliton oscillation frequency 
on the system parameters can be very useful for experimental studies of relative 
contribution of multi-particle collisions into BEC dynamics by varying the value of g by 
means of Feshbach resonance. Moreover, it gives the possibility to measure background 
densities or scattering lengths directly from the dark soliton dynamics. The possibility 
to manage soliton oscillations via temporal changes of the scattering lengths is also 
an attractive feature emerging from our study which could be of interest for practical 
applications. 

The paper is organized as follows. In Section II we discuss the general theory of 
the dynamics of a dark soliton of the NLS equation with an arbitrary nonlinearity. 
The results will be used to study the cases of both deep and shallow soliton dynamics. 
In Section III we consider oscillations of dark solitons for the specific case of a cubic- 
quintic nonlinearity. We consider the cases of a repulsive cubic nonlinearity and either 
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repulsive or attractive quintic nonlinearity, the last being done for parameters satisfying 
the condition of modulation stability of BEC. In Section IV analytical results will be 
compared with direct numerical simulations of the cubic-quintic GP equation for the 
case of repulsive and attractive cubic and quintic nonlinearity, respectively. In the last 
Section we draw our conclusions and summarize the main results of the paper. 



2. Motion of a dark soliton in a trap: general theory 

First, we have to find the soliton solution of the generalized NLS equation (jlj) with zero 
potential V(x) = but under the boundary condition that the condensate density p(x) 
tends to the limiting background density po as |x| — > oo and the condensate's velocity 
vanishes at infinity. To this end, it is convenient to look for the solution in the form 



V> = v^expMO-^O (11) 

where 

» = f(po) (12) 

is the chemical potential of the condensate and its density p and the local flow velocity 
u = (f x depend only on the variable £ = x — vt with v being the soliton velocity. The 
formulated above boundary conditions take the form 

P -> Po, n °> as If I 00 . (13) 
Substitution of (ITT]) into (j3J) with V(x) = and separation of real and imaginary parts 
yields the system of two equations for and = <p%. One of these equations can 
be integrated at once to give 

u(0 = tpt = v(l - ^\ (14) 

where the integration constant is chosen according to the conditions ( fl3l) . Then the 
second equation can be transformed to 

IpI - \pp& + p 2 f(p) - + \ApI - p 2 ) = o (15) 

and this equation can also be integrated once. As a result we obtain the equation 

P\ = Q(P) (16) 

where 

rpo 

Q(p) = 8p [p- f(p')]dp' - 4v 2 (p - pf (17) 
J P 

and again the integration constant is chosen according to the condition (fl~3|) . Thus, we 
have reduced the problem of finding the density profile to the inversion of the integral 

^ Lm y/Q(p) ^ ^ 

where p m is the minimal density at the center £ = of the soliton. The function 
Q(p) has a double zero at p = p and the condition dQ(p)/dp\ = reproduces the 
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relation (|T2"|) . At p = p m the function Q(p) has a single zero, and this condition yields 
the relationship between p m and the soliton velocity v, 

2 Qo(Pm) 



4(p - Pm) 2 ' 



where 



Qo(p) = 8p / [/(p ) - f{p')W 



(19) 



(20) 



is the function (JTTj) with zero value of the velocity w and we have used here Eq. ffT2l for 
p. 

Now we have to obtain the expression for the energy of the dark soliton. To this 
end, we notice that the function ^ = yj~p exp(itp) = ip ■ exp(ipt) satisfies the equation 

iU** + (/W(l*| 2 ))* = (21) 



which can be written in a Hamiltonian form 
SH 

i% - 



with 



H 



ii*. 



+ 



l*l 2 



f{p)dp-p\m\ 



dx. 



(22) 



(23) 



To get the contribution of the dark soliton into the energy of BEC, we subtract the 
background energy H = [J^° f(p)dp — p\^\ 2 ] dx and express the resulting formula 
in terms of p and tp x , 



E = H — Ho 



2 I Px 



p 



[/(Po) " f{ P ')W 



dx. (24) 



Substitution of expressions (fT4|) and (|T6l) for the dark soliton solution casts this formula 
into 

fOO fpo 

[/ (Po) - f(p'Wdx. (25) 



OO J p 



At last, integration with respect to x can be transformed to integration with respect to 
p from p rn to p with account of (fT8|) . 



E 



r»° dp r [f{po) _ f{p , )]dp , = i r° go(p) ip 

A™ y/Q(p) Jo 2 A™ P\/Q(p) 



pVQipj 



(26) 



Here p m can be considered according to fTT9~j) as a function of p an d hence the soliton 
energy is also a function of these two variables, 

E = E(p ,v 2 ). (27) 

Now we turn to the problem of motion of the dark soliton along a nonuniform 
condensate confined in a trap with the axial potential V(x). Then in the Thomas- Fermi 
(TF) approximation the density p(x) is a function of x determined implicitly by the 
equation 

p-V(x) = f(p(x)). (28) 
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When this TF distribution of the density p = p(x) is substituted into (127)) instead of p , 
we obtain the energy of the soliton moving with velocity v and located at this moment 
at the point x. As was shown in Ref. [15], the velocity v changes during the motion of 
the dark soliton in such a way that the soliton energy (I27j) is conserved. If we denote 
the soliton coordinate as x — X(t), then v = X(t) and Eq. (1271) converts into 

E(p(X),X 2 ) = E , (29) 

where E is the initial energy of the soliton. Since the function E(p, v 2 ) is known, 
Eq. ( 1291) is the differential equation for finding the soliton coordinate X as a function 
of time t. In general, this first order ordinary differential equation should be solved 
numerically. Here we shall consider the most interesting limiting cases admitting 
complete analytical study. 



2.1. Small amplitude oscillations of a deep soliton 
Let us consider oscillations of a deep soliton, 

Pm < PO, (30) 

where po is the maximal density of the condensate at its "top" when it is confined in 
a trap with the potential V(x) (see ()28l) ). In this case the amplitude of oscillations is 
much less than the TF radius of the condensate and the velocity of its motion is much 
less than the local sound velocity 

v < c s , c s = y/pof'(po)- (31) 

Therefore, if we expand the energy into series with respect to small amplitude and 
velocity, then the first nontrivial terms in ()29l) will give the expression 

^m,X 2 + V(X) = const, (32) 

where V(X) changes little for the small amplitude vibrations. Hence, we get the "energy 
conservation law" for a Newtonian particle moving in the potential V(X), where m* is 
the effective mass of the dark soliton. 
For small v we obtain from (fT9"j) 

Pm = ^ (33) 

2a 



where 



Po 

[/(Po) - f(p)]dp. (34) 



o 



To estimate the first term of the series expansion of E with respect to v 2 , we split the 
integral ()26j) into two terms (see [T6]). 

E — Ei + E2, 

E = 1 rPi _Q^ d E = 1 p _Q^ d (35) 
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where we have introduced an intermediate integration limit pi so that 

p m < pi < p - (36) 

Dependence of E\ on v 2 results mainly from the lower integration limit ( 133]) and in the 
main approximation we get 



E 1 



pi 



-dp 



V0(P)-- (37) 

' Pm P 

Now, both limits of integration are small and p m is a simple zero of the function Q(p), 
hence we can approximate it as 

dQ(p) 



Q(p) 



dp 



(p- Pm) = 8d(p-p m ), 



(38 



where a is defined in ( 13"4"1) . Then substitution of ( 1381) into (1371) and elementary integration 
yield 

Ei = v / 8a{p 1 



Pr 



(39) 



We need the coefficient of v in the expansion of E in powers of v and the contribution 
of Ei into this coefficient is equal to 

dEhdprn = Pq , 4Q , 
dp m dv 2 \J2api 

As we see, it diverges in the limit pi — > 0. To cancel this divergent term, we have to 
identify a similar term in the derivative dE 2 /dv 2 given by the expression 



dEo 



dv 1 



v=Q 



Pi 



(Po - P? 

pVQo(p) 



dp. 



(41) 



According to Eq. (140jl . the diverging term is proportional to p 1 . Correspondingly, we 



transform the integral in (14ip in the following way by integration by parts, 



(42) 



dE 2 _ f 

dv 2 \v=a J i 



Po {po-p) 2 



P 1 P 3/2 y/Qa(p)/p 
2(P0-P) 2 



y/Qo(fi) 



dp. 



rpo 2 d 
J pi ^fp dp 



Pi 



y/Qo(fi) 



The first term vanishes at the upper limit p = po and at the lower limit p = pi <C po 
we can approximate Qo(p) as Qo(p) — 8ap (see ( )38|) ). In the second term the integral is 
convergent and we can take here the lower integration limit equal to pi = 0. As a result 

we get 

Pi + 2 f P0 [ ^(po-p) 2 " 
v/SopT Jo y/pdp [ ^Qo(p) 

The first term cancels with the contribution ( j40l) and we obtain the expression for the 
soliton mass 



dEo 



dv 2 



v=Q 



(43) 



dE 7 



dv 2 



v=0 



rp0 dp d 

yfpdp 



yp(po - pf 
Vo^pj 



(44) 
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for the motion of a dark soliton near the "top" of the condensate. 

Introduction of the potential V(x) leads to the shift of the chemical potential 
fi — > fi — V(x). Taking into account the thermodynamic formula 

* = -f <«> 

and the expression for the energy in the limit v 2 — > 0, 

dp 

ff ' ,i/ 2 (46) 

^- 2 f (47) 
Jo Vvo(p) 

we find the expansion of the soliton energy as 



e = e ° =h w -7s=/;>-/(<>w 



hence 



E = E Q + 



dE 



■v 2 + ^- (-V(x)) = E + \m s v 2 + NV(x), (48) 
v 2 =0 an z 



dv 2 

and division by N yields the energy conservation law for the soliton's motion 



where 



— V + V(x) = const (49) 

m * = jj L ( 50 ) 

is the effective mass of the soliton's motion. Hence, in case of a harmonic trap (jSJ) we 
get the formula 

uj = ^ (51) 

nru 



for the frequency of oscillations of deep solitons. 
For a usual GP equation we have 

f{p) = P, (52) 

hence 

Qo(p) = 4p(p - p)\ m s = -Ay/ft, N = -2^ (53) 

and 

m* = 2 (54) 
so that (l5"Tj) reproduces the well-known result. 



Dark soliton oscillations 



10 



2.2. Motion of a shallow soliton 

In the shallow soliton limit 

Po - P < Po 
the series expansion of Eq. (fl9l) gives 

_ ,2 rsj 



Pi 



or with account of the expression ( TlQl) for the sound velocity we obtain 



Pr, 



Po 



(55) 



(56) 



(57) 



/'(po) + fpo/"(po)' 

In this limit the dark soliton can be considered in the Korteweg-de Vries (KdV) 
approximation (see, e.g. [21]) for which the NLS equation with V(x) = reduces to the 
equation 



P't + C sp'x + 



2c e 



/'(Po) + ^Po/"(po) 



P Pz - Pxxx = 

OCe 



(58 



for evolution of small deviations p' = po — p of the density from the background value 
Po- This equation has the soliton solution 



1 



P = Po - P 



/'(Po) + |po/"(po) cosh 2 [ K (x-vt)} 



where 



K 



c s (c 2 s - V 2 ) 



(59) 



(60) 



/'(Po) + |po/"(po) 

is the inverse soliton's half-width and the minimal density p m at the soliton's center 
x = vt is given by the expression ([3] 
In the same approximation 



Q{p) = 4 



/'(Po) + 3P0/" (Po) 



(P0-P) 2 (p-Pm) 



and calculation of the soliton energy (1261) gives 
4 f'(po) 



E 



3 [/'(po) + |po/"(po 



2 _ t; 2 )3 /2_ 



(61) 



(62) 



If the width 1/k of the soliton is much less than the condensate's size (the TF 
radius), then we can consider the motion of such a soliton in the trap by replacement 
of po by the solution p(x) of Eq. ( 1281) and of v by its variable velocity dX/dt to get the 
equation 



dXV , J 3E 

^ p/(p) - \mp) 



Z'(P) + g/"(p) 



2/3 



(63) 
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where p = p(X) is given by Eq. (j28p . For example, in case of / = hp 2 , h > 0, we have 
f'(p) = 2hp and p — V(X) = hp 2 , hence f )63l) reduces to 

As we see, for a harmonic trap (jSJ) the frequency of oscillations depends on the energy 
E of the soliton. However, if E is small enough, 

E « (65) 

then the last term in the left-hand size can be neglected and we return to the usual energy 
conservation law for a Newtonian particle which for a harmonic trap gives Eq. (jSJ) in 
agreement with [TH | 



3. Oscillation of a dark soliton for the cubic-quintic nonlinearity 

Now we shall consider in some detail the motion of a dark soliton described by the 
generalized NLS equation with /(p) given by Eq. (jSJ), i.e., 

f(p)=gp + hp 2 . (66) 

Conditions for realization of motion of a shallow soliton are very restrictive and we shall 
confine ourselves to the most practical case of small oscillations of a deep soliton (Sec. 
2A). 

Substitution of ( 1661) into ( 1201) gives 



Qo(p) = 8p(p - p) 



-g + hi{2 Po + p) 



(67) 



Situations with h > and h < should be considered separately. 
3.1. Repulsive quintic interaction (h > 0) 

In this case elementary integration of Eq. (fl8|) with Q(p) given by (|67|) yields the profile 
of the density in the dark solitons solution 

ic\ Qo - p+)(po - p-) , . 

p{0 = Po-~, : 7T ( 77 , ^ = x-ut, (68) 

(P+ - P-) cosh «) + p - p+ 

where 



and 



« = V / 2MPo-p+)(po-P-)/3. (70) 
Substitution of (167|) into (147|) and simple calculation give 



N- .Tin ( / 2, 'P" | % + 2/W, 
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X t 

Figure 1. Left panel. Time evolution of a deep dark soliton in a BEC with repulsive 
two- and three-body interactions obtained from numerical simulation of the cubic- 
quintic GPE with g = 1.0, h = 0.5 and parabolic potential V(x) = \uq with uiq — 0.894 
(the potential is depicted at the top of the figure for scale comparison). Right 
panel. Dark soliton center coordinate (bottom curve) and soliton density background 
corresponding to the dynamics depicted in the right panel. The numerical oscillation 
frequency of the soliton is uj = 0.607, while the time averaged background density, 
indicated by the dashed line in the figure, is (po) = 3.674. The thin dotted line at the 
bottom refers to the dynamics of the center of mass of the whole condensate. 



Similar calculation of (PHI) yields 

2hp l3(g + 2hp )\ ^6hp (g + 2hp ) 



'6 
i/ - 



In 



3g + 4hp y 3g + 4hp ) 3g + 4hp 

Hence the effective mass of the soliton motion is equal to 



.(72) 



m. = 1 + y/Mpoig + Zhpo) (73) 

(3^ + ^ )ln(v^+yigg) 

For hp g this formula reproduces 

m* = 2, hp < g, (74) 

in agreement with (jSJ), and for hp ^> g we get 

m * = 1 H TTT^' h Po > ( 75 ) 

21ni±^ 

in agreement with (jUj). 

To check these results, we have numerically integrated the GP equation for this 
case using as initial condition the exact dark soliton solution of the cubic-quintic NLS 
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Figure 2. Dependence of the ratio luq/uj = y/m* on the BEC background density po for 
repulsive interactions, g — 1.0, h = 0.5, (left panel) and for opposite sign interactions 
g = 1.0, h = —0.5 (right panel). Continuous curves refer to the analytical results 
in Eqs. ([751) . Q82[) . while the dots are obtained from direct numerical solution of the 
generalized NLS equation (TJ} with f(p) given by Eq. (l6Gl) . 



in Eq. ([T]). The parabolic potential was introduced adiabatically by rising with time 
the frequency uj from zero (pure cubic-quintic solution) up to a desired value. Effects 
of adiabatic changes of parameters on BEC bright and dark solitons were considered in 
Ref. [22] and the effectiveness of the method was demonstrated. In Fig. [1] we depict the 
dynamics of a repulsive condensate with its dark soliton inside in presence of a parabolic 
trap raised adiabatically as described above. In the right panel of this figure we show the 
dynamics of the center position X(t) of the dark soliton and of the center of mass of the 
whole condensate X cm = J x\^\ 2 dx together with the time evolution of the background 
density at the soliton position. We see that both the position of the soliton and of 
the BEC center of mass are very regular and relatively small amplitude excitations 
develop on the background during the evolution. Note that the center of mass of the 
condensate oscillates with the trap frequency in agreement with Ehrenfest theorem, 
while the oscillation frequency and the effective mass of the dark soliton are in very good 
agreement with our analytical predictions u = 0.599, m* = 2.227 (for the comparison 
with numerics we take as po the time average of the soliton background measured during 
the evolution). A more extended comparison between the theory and numerical results 
is presented in the left panel of Fig. [2] in which we see that in general there is a good 
qualitative agreement which becomes also quantitatively good at higher values of soliton 
backgrounds. The deviation observed at lower values of po can be ascribed, probably, to 
a large soliton width comparable with the size of the whole condensate. In this situation 
the conditions of validity of the theory are not satisfied well enough. 
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3.2. Attractive quintic interaction (h < 0) 

For h < the integrand function in Eq. (1471) is real for 



" < 4' (76) 

and this inequality coincides with the condition of modulation stability of BEC, namely 
the condition that the sound velocity c s = a/ Pof'(po) is rea l. It is worth noticing that 
this condition is more restrictive than the condition p < d/\h\ that the nonlinearity 
function is positive, /(p) = gp — \h\p 2 > 0. Obviously, we can discuss the soliton 
oscillations for BECs satisfying the condition (1761) only. 

In this case Q(p) has the double zero po and the simple zeroes p + > po and 
< p_ < p where 




Hence there exist both dark and bright soliton solutions. We shall consider here motion 
only of the dark soliton with the density profile 

(po ~ P~)(p+ ~ Po) 
(p+ ~ P-) cosh 2 «) - (p - p_ 
where 

k = v / 2|/ i |(p -p_)(p + -po)/3. (79) 
Now calculation of the number of particles N and the soliton mass m s gives 



tc\ kho - h- \y+ ho t TQ\ 

P(0 = PO- — \ 2/ -y, ~ — ' ( 78 ) 




iV = -J^arcsinJ o ^ (80) 



\h\ y3g-A\h\p 



and 



m. 




arcsm 



2\h\p , 2^6\h\ Po {g-2\h\ Po ) 



3g-A\h\p 3g-A\h\po 
Hence the effective mass of the soliton motion is equal to 



(81) 



m, = I2i = l + VQ\h\Po(9-2\h\Po) (g2) 

N (3,-4|/ i |p )arcsin v ^^: 

The condition ( 1761) means that in these expressions the argument of the arcsin function 
is less than unity. (Of course, expressions (l80l - (!82l) can be obtained from (17Tl) - (!73l by 
means of their analytical continuation to negative values of h.) 

For \h\p <C g we again reproduce (174]) . as it should be. For g — 2\h\p <C g we find 
the limiting behavior 

2 

m*=H —\Jg- 2|/i|p , p - 2\h\p <C p. (83) 

Thus, function of po has a branching point singularity at po c = g/(2|/i|) where 

the system becomes modulationally unstable. 
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The plot of the ratio of the trap frequency to the soliton's oscillations frequency is 
illustrated in Fig. 2 (see right panel) and it demonstrates a very good agreement with 
the results of the numerical simulations. 

4. Conclusion 

In this paper we have considered the dynamics of a dark soliton inside a BEC described 
by the quasi-one dimensional GP equation with a parabolic trapping potential. In 
particular, we have developed the theory which permits one to study the motion of a 
dark soliton through the stationary background state of BEC and to derive the formulae 
for the frequency of the oscillation of deep solitons which are very effective for arbitrary 
form of the nonlinearity function f(p). When the form of f(p) is appropriately specified, 
the formulae for the frequency of oscillations reproduce the results of all previously 
studied special cases. We applied our theory to the study of the practically interesting 
case of cubic-quintic nonlinearity function f(p) = gp + hp 2 and we have demonstrated, 
as a new effect, the dependence of the frequency of the oscillations on the background 
density of BEC. In principle this effect can be used for the experimental study of higher 
nonlinearity on BEC's dynamics. 
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